Deep dive regional communities - using relative abundance prox

source('../env.R')

Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
── Attaching core tidyverse packages ────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.4.4     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     ── Conflicts ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errorsRegistered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      

Attaching package: ‘dbplyr’

The following objects are masked from ‘package:dplyr’:

    ident, sql

Loading required package: ape

Attaching package: ‘ape’

The following object is masked from ‘package:dplyr’:

    where

Loading required package: maps

Attaching package: ‘maps’

The following object is masked from ‘package:purrr’:

    map

Please cite the eBird Status and Trends data using: 
  Fink, D., T. Auer, A. Johnston, M. Strimas-Mackey, S. Ligocki, O. Robinson, 
  W. Hochachka, L. Jaromczyk, C. Crowley, K. Dunham, A. Stillman, I. Davies, 
  A. Rodewald, V. Ruiz-Gutierrez, C. Wood. 2023.
  eBird Status and Trends, Data Version: 2022; Released: 2023. Cornell Lab of
  Ornithology, Ithaca, New York. https://doi.org/10.2173/ebirdst.2022

This version of the package provides access to the 2022 version of the eBird
Status and Trends Data Products. Access to the 2022 data will be provided 
until November 2024 when it will be replaced by the 2023 data. At that 
point, you will be required to update this R package and transition to using 
the new data.
terra 1.7.71

Attaching package: ‘terra’

The following object is masked from ‘package:phytools’:

    rescale

The following objects are masked from ‘package:ape’:

    rotate, trans, zoom

The following object is masked from ‘package:tidyr’:

    extract

Loading required package: permute
Loading required package: lattice
This is vegan 2.6-4

Attaching package: ‘vegan’

The following object is masked from ‘package:phytools’:

    scores

Loading required package: mvtnorm
Loading required package: survival
Loading required package: TH.data
Loading required package: MASS

Attaching package: ‘MASS’

The following object is masked from ‘package:terra’:

    area

The following object is masked from ‘package:dplyr’:

    select


Attaching package: ‘TH.data’

The following object is masked from ‘package:MASS’:

    geyser


Attaching package: ‘ggpubr’

The following object is masked from ‘package:terra’:

    rotate

The following object is masked from ‘package:ape’:

    rotate

Loading required package: nlme

Attaching package: ‘nlme’

The following object is masked from ‘package:dplyr’:

    collapse

Registered S3 method overwritten by 'GGally':
  method from   
  +.gg   ggplot2

Attaching package: ‘GGally’

The following object is masked from ‘package:terra’:

    wrap

ggtree v3.10.0 For help: https://yulab-smu.top/treedata-book/

If you use the ggtree package suite in published research, please cite the appropriate paper(s):

Guangchuang Yu, David Smith, Huachen Zhu, Yi Guan, Tommy Tsan-Yuk Lam. ggtree: an R package for visualization and annotation of phylogenetic
trees with their covariates and other associated data. Methods in Ecology and Evolution. 2017, 8(1):28-36. doi:10.1111/2041-210X.12628

Guangchuang Yu. Using ggtree to visualize data on tree-like structures. Current Protocols in Bioinformatics. 2020, 69:e96. doi:10.1002/cpbi.96

Guangchuang Yu.  Data Integration, Manipulation and Visualization of Phylogenetic Trees (1st edition). Chapman and Hall/CRC. 2022,
doi:10.1201/9781003279242 

Attaching package: ‘ggtree’

The following object is masked from ‘package:nlme’:

    collapse

The following object is masked from ‘package:ggpubr’:

    rotate

The following objects are masked from ‘package:terra’:

    flip, inset, rotate

The following object is masked from ‘package:ape’:

    rotate

The following object is masked from ‘package:tidyr’:

    expand


Attaching package: ‘foreach’

The following objects are masked from ‘package:purrr’:

    accumulate, when


Attaching package: ‘scales’

The following object is masked from ‘package:terra’:

    rescale

The following object is masked from ‘package:phytools’:

    rescale

The following object is masked from ‘package:purrr’:

    discard

The following object is masked from ‘package:readr’:

    col_factor

Spherical geometry (s2) switched off

Species in communities

It seems reasonable to expect that cities with simialr regional pools will have similar species entering the city, and thus a similar response to urbanisation.

Load data

city_effort = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'city_effort.csv'))

Rows: 342 Columns: 7── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): city_name
dbl (6): city_id, total_city_checklists, total_city_locations, total_city_effort, total_city_area_m2, percentage_total_city_area_surveyed
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

city_effort

communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))

Rows: 2462 Columns: 12── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, jetz_species_name, seasonal, presence, origin
dbl (4): city_id, distance_to_northern_edge_km, distance_to_southern_edge_km, relative_abundance_proxy
lgl (3): present_urban_high, present_urban_med, present_urban_low
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

communities_summary = communities %>% group_by(city_id) %>% summarise(
  regional_pool_size = n(), 
  urban_pool_size = sum(relative_abundance_proxy > 0)
) %>% left_join(city_effort %>% dplyr::select(city_id, percentage_total_city_area_surveyed))

Joining with `by = join_by(city_id)`

ggplot(communities %>% filter(relative_abundance_proxy > 0), aes(x = relative_abundance_proxy)) + geom_bar(stat = "bin")

city_points = st_centroid(read_sf(filename(CITY_DATA_OUTPUT_DIR, 'city_selection.shp')))

Warning: st_centroid assumes attributes are constant over geometriesWarning: st_centroid does not give correct centroids for longitude/latitude data

community_data_metrics = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'community_assembly_metrics_using_relative_abundance.csv')) %>%
  dplyr::select(city_id, mntd_normalised, fdiv_normalised, mass_fdiv_normalised, kipps_distance_fdiv_normalised, trophic_trait_fdiv_normalised, gape_width_fdiv_normalised) %>%
  left_join(read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv'))) %>%
  left_join(communities_summary) %>%
  left_join(city_points[,c('city_id', 'city_nm')]) %>%
  rename(city_name='city_nm') %>%
  na.omit() %>%
  arrange(city_id)

Rows: 341 Columns: 43── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (43): mntd_normalised, mntd_actual, mntd_min, mntd_max, mntd_mean, mntd_sd, fdiv_normalised, fdiv_actual, fdiv_min, fdiv_max, fdiv_mean, fdiv_sd, mass_fdi...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Rows: 342 Columns: 2── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): core_realm
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`

community_data_metrics

test_value = function(name, normalised_list) {
  wilcox_test_result = wilcox.test(normalised_list, mu = 0.5)
  
  significance = ifelse(wilcox_test_result$p.value < 0.0001, '***', 
                        ifelse(wilcox_test_result$p.value < 0.001, '**', 
                               ifelse(wilcox_test_result$p.value < 0.01, '*', '')))
  m = mean(normalised_list)
  
  paste(name, 'mean', round(m, 2), significance)
}

test_value('Global MNTD', community_data_metrics$mntd_normalised)

[1] "Global MNTD mean 0.5 "

test_value('Global beak gape', community_data_metrics$gape_width_fdiv_normalised)

[1] "Global beak gape mean 0.59 ***"

test_value('Global kipps distance', community_data_metrics$kipps_distance_fdiv_normalised)

[1] "Global kipps distance mean 0.75 ***"

Load trait data

traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_jetz.csv'))

Rows: 304 Columns: 6── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): jetz_species_name
dbl (5): gape_width, trophic_trait, locomotory_trait, mass, kipps_distance
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(traits)

fetch_normalised_traits = function(required_species_list) {
  required_traits = traits %>% filter(jetz_species_name %in% required_species_list)
  
  required_traits$gape_width_normalised = normalise(required_traits$gape_width, min(required_traits$gape_width), max(required_traits$gape_width))
  required_traits$kipps_distance_fdiv_normalised = normalise(required_traits$kipps_distance, min(required_traits$kipps_distance), max(required_traits$kipps_distance))
  required_traits$mass_normalised = normalise(required_traits$mass, min(required_traits$mass), max(required_traits$mass))
  
  traits_normalised_long = required_traits %>% pivot_longer(cols = c('gape_width_normalised', 'kipps_distance_fdiv_normalised', 'mass_normalised'), names_to = 'trait', values_to = 'normalised_value') %>% dplyr::select(jetz_species_name, trait, normalised_value)
  traits_normalised_long$trait = factor(traits_normalised_long$trait, levels = c('gape_width_normalised', 'kipps_distance_fdiv_normalised', 'mass_normalised'), labels = c('Gape Width', 'Kipps Distance', 'Mass'))
  
  traits_normalised_long
}

fetch_normalised_traits(c('Aplopelia_larvata', 'Chalcophaps_indica', 'Caloenas_nicobarica'))

Read in our phylogeny

phylo_tree = read.tree(filename(TAXONOMY_OUTPUT_DIR, 'phylogeny.tre'))
ggtree(phylo_tree, layout='circular')

Load resolve ecoregions

resolve = read_resolve()

Warning: st_buffer does not correctly buffer longitude/latitude datadist is assumed to be in decimal degrees (arc_degrees).
Warning: st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees

Create helper functions

to_species_matrix = function(filtered_communities) {
  filtered_communities %>% 
    dplyr::select(city_id, jetz_species_name) %>% 
    distinct() %>%
    mutate(present = TRUE) %>% 
    pivot_wider(
      names_from = jetz_species_name, 
      values_from = "present", 
      values_fill = list(present = F)
    ) %>% 
    tibble::column_to_rownames(var='city_id')
}

community_nmds = function(filtered_communities) {
  species_matrix = to_species_matrix(filtered_communities)
  nmds <- metaMDS(species_matrix, k=2, trymax = 30)
  nmds_result = data.frame(vegan::scores(nmds)$sites)
  nmds_result$city_id = as.double(rownames(nmds_result))
  rownames(nmds_result) = NULL
  nmds_result
}

https://www.datacamp.com/tutorial/k-means-clustering-r

scree_plot = function(community_nmds_data) {
  # Decide how many clusters to look at
  n_clusters <- min(10, nrow(community_nmds_data) - 1)
  
  # Initialize total within sum of squares error: wss
  wss <- numeric(n_clusters)
  
  set.seed(123)
  
  # Look over 1 to n possible clusters
  for (i in 1:n_clusters) {
    # Fit the model: km.out
    km.out <- kmeans(community_nmds_data[,c('NMDS1','NMDS2')], centers = i, nstart = 20)
    # Save the within cluster sum of squares
    wss[i] <- km.out$tot.withinss
  }
  
  # Produce a scree plot
  wss_df <- tibble(clusters = 1:n_clusters, wss = wss)
   
  scree_plot <- ggplot(wss_df, aes(x = clusters, y = wss, group = 1)) +
      geom_point(size = 4) +
      geom_line() +
      geom_hline(linetype="dashed", color = "orange", yintercept = wss) +
      scale_x_continuous(breaks = c(2, 4, 6, 8, 10)) +
      xlab('Number of clusters')
  scree_plot
}

cluster_cities = function(city_nmds, cities_community_data, centers) {
  set.seed(123)
  kmeans_clusters <- kmeans(city_nmds[,c('NMDS1', 'NMDS2')], centers = centers, nstart = 20)
  city_nmds$cluster = kmeans_clusters$cluster
  cities_community_data %>% left_join(city_nmds) %>% mutate(cluster = as.factor(cluster))
}

plot_nmds_clusters = function(cluster_cities) {
  cluster_cities %>% dplyr::select(city_id, city_name, NMDS1, NMDS2, cluster) %>% distinct() %>%
  ggplot(aes(x = NMDS1, y = NMDS2, colour = cluster)) + geom_point() + geom_label_repel(aes(label = city_name))
}

get_presence_cell_width = function(city_cluster_data_metrics) {
  10 * length(unique(city_cluster_data_metrics$city_id))
}

get_presence_cell_height = function(city_cluster_data_metrics) {
  species = species_in_cluster = communities %>% 
    filter(city_id %in% city_cluster_data_metrics$city_id) %>% 
    dplyr::select(jetz_species_name) %>% 
    distinct()
  
  10 * nrow(species)
}

city_metric_height = 30
traits_width = 50
phylo_tree_width = 125
title_height = 8

get_image_height = function(city_cluster_data_metrics) {
  get_presence_cell_height(city_cluster_data_metrics) + (3 * city_metric_height) + title_height
}

get_image_width = function(city_cluster_data_metrics) {
  get_presence_cell_width(city_cluster_data_metrics) + traits_width + phylo_tree_width
}

species_in_city_label = function(species) {
  paste(
    ifelse(species$seasonal == 'Resident', '', substring(species$seasonal, 0, 1)),
    ifelse(species$origin == 'Native', '', substring(species$origin, 0, 1)),
    ifelse(species$distance_to_northern_edge_km > 200, '', paste('NRL', round(species$distance_to_northern_edge_km), sep = ' ')),
    ifelse(species$distance_to_southern_edge_km > 200, '', paste('SRL', round(species$distance_to_northern_edge_km), sep = ' ')),
    sep = '\n'
  )
}

species_in_city_label(head(communities))

[1] "\nI\n\n" "\nI\n\n" "\nI\n\n" "P\n\n\n" "\n\n\n"  "\nI\n\n"

plot_city_cluster = function(city_cluster_data_metrics, title) {
  species_in_cluster = communities %>% 
    filter(city_id %in% city_cluster_data_metrics$city_id) %>% 
    dplyr::select(jetz_species_name, city_name, relative_abundance_proxy, seasonal, origin, distance_to_northern_edge_km, distance_to_southern_edge_km)
  
  species_in_cluster$label = species_in_city_label(species_in_cluster)
  tree_cropped <- ladderize(drop.tip(phylo_tree, setdiff(phylo_tree$tip.label, species_in_cluster$jetz_species_name)))
    
  gg_tree = ggtree(tree_cropped)
  
  gg_presence = ggplot(species_in_cluster, aes(x=city_name, y=jetz_species_name)) + 
          geom_tile(aes(fill=relative_abundance_proxy)) + 
          geom_text(aes(label=label), size=0.75) +
          scale_fill_gradientn(colours=c("#98FB98", "#FFFFE0", "yellow", "orange", "#FF4500", "red", "red"), values=c(0, 0.00000000001, 0.1, 0.25, 0.5, 0.75, 1), na.value = "transparent") +
          theme_minimal() + xlab(NULL) + ylab(NULL) + 
          theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) + 
          labs(fill='Urban Proxy Abundance')
  
  species_in_cluster_traits = fetch_normalised_traits(species_in_cluster$jetz_species_name)
  
  gg_traits = ggplot(species_in_cluster_traits, aes(x = trait, y = jetz_species_name, size = normalised_value)) + geom_point() + theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text.y=element_blank()) + xlab(NULL) + ylab(NULL) + labs(size = "Normalised Value")
  
  gg_cities_mntd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = mntd_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("MNTD") + ylim(0, 1)
  
  gg_cities_gape_fd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = gape_width_fdiv_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("Gape") + ylim(0, 1)
  
  gg_cities_loco_fd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = kipps_distance_fdiv_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("Kipps") + ylim(0, 1)
  
  gg_title = ggplot() + labs(title = title, subtitle = paste(
    test_value('MNTD', city_cluster_data_metrics$mntd_normalised),
    test_value('FDiv', city_cluster_data_metrics$fdiv_normalised),
    test_value('Kipps distance', city_cluster_data_metrics$kipps_distance_fdiv_normalised),
    test_value('Gape width', city_cluster_data_metrics$gape_width_fdiv_normalised),
    sep = '\n'
  )) + theme_minimal() + theme(plot.subtitle=element_text(size=8, hjust=0, color="#444444"))
  
  gg_presence_height = get_presence_cell_height(city_cluster_data_metrics)
  gg_presence_width = get_presence_cell_width(city_cluster_data_metrics)
  
  gg_presence %>% insert_top(gg_cities_loco_fd, height = (city_metric_height / gg_presence_height)) %>% insert_top(gg_cities_gape_fd, height = (city_metric_height / gg_presence_height)) %>% insert_top(gg_cities_mntd, height = (city_metric_height / gg_presence_height)) %>% insert_left(gg_tree, width = (phylo_tree_width / gg_presence_width)) %>% insert_right(gg_traits, width = (traits_width / gg_presence_width)) %>% insert_top(gg_title, height = (title_height / gg_presence_height))
}

REGION_DEEP_DIVE_FIGURES_OUTPUT = mkdir(FIGURES_OUTPUT_DIR, 'appendix_regional_deep_dive_using_abundance')

Nearctic

nearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Nearctic')
nearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
 [1] "San Jose"                 "Los Angeles"              "Concord"                  "Tijuana"                  "Bakersfield"             
 [6] "Fresno"                   "Sacramento"               "Mexicali"                 "Hermosillo"               "Las Vegas"               
[11] "Phoenix"                  "Tucson"                   "Durango"                  "Portland"                 "Chihuahua"               
[16] "Aguascalientes"           "Seattle"                  "Ciudad Juárez"            "San Luis Potosí"          "Mexico City"             
[21] "Saltillo"                 "Vancouver"                "Salt Lake City"           "Albuquerque"              "Monterrey"               
[26] "Nuevo Laredo"             "San Antonio"              "Denver"                   "Austin"                   "Houston"                 
[31] "Dallas"                   "Oklahoma City"            "Calgary"                  "New Orleans"              "Kansas City"             
[36] "Omaha"                    "St. Louis"                "Bradenton"                "Tampa"                    "Minneapolis [Saint Paul]"
[41] "Atlanta"                  "Orlando"                  "Louisville"               "Chicago"                  "Indianapolis"            
[46] "Milwaukee"               

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

nearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% nearctic_cities_community_data$city_id)) 

Run 0 stress 0.1005678 
Run 1 stress 0.1017503 
Run 2 stress 0.1000046 
... New best solution
... Procrustes: rmse 0.007231564  max resid 0.03470212 
Run 3 stress 0.1000046 
... Procrustes: rmse 0.000004044592  max resid 0.00001391239 
... Similar to previous best
Run 4 stress 0.1205431 
Run 5 stress 0.1012776 
Run 6 stress 0.122422 
Run 7 stress 0.1005678 
Run 8 stress 0.1000046 
... Procrustes: rmse 0.000006352527  max resid 0.00001951675 
... Similar to previous best
Run 9 stress 0.1217145 
Run 10 stress 0.1000046 
... New best solution
... Procrustes: rmse 0.000003781293  max resid 0.00001247782 
... Similar to previous best
Run 11 stress 0.1217437 
Run 12 stress 0.1217437 
Run 13 stress 0.1000046 
... Procrustes: rmse 0.000004334204  max resid 0.00001530842 
... Similar to previous best
Run 14 stress 0.1000046 
... Procrustes: rmse 0.000004734417  max resid 0.00001660964 
... Similar to previous best
Run 15 stress 0.1012776 
Run 16 stress 0.1235457 
Run 17 stress 0.1000046 
... Procrustes: rmse 0.000004211501  max resid 0.00001320069 
... Similar to previous best
Run 18 stress 0.1005678 
Run 19 stress 0.1005678 
Run 20 stress 0.1017503 
*** Best solution repeated 4 times

nearctic_cities_nmds

scree_plot(nearctic_cities_nmds)

nearctic_cities = cluster_cities(city_nmds = nearctic_cities_nmds, cities_community_data = nearctic_cities_community_data, centers = 4)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(nearctic_cities)

nearctic_biomes = st_crop(resolve[resolve$REALM == 'Nearctic',c('REALM')], xmin = -220, ymin = 0, xmax = 0, ymax = 70)

although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning: attribute variables are assumed to be spatially constant throughout all geometries

 
ggplot() + 
  geom_sf(data = nearctic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = nearctic_cities, aes(geometry = geometry, color = cluster))

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_clusters.jpg'))

Saving 7.29 x 4.51 in image

test_value('Nearctic MNTD', nearctic_cities$mntd_normalised)

Warning: cannot compute exact p-value with ties

[1] "Nearctic MNTD mean 0.6 "

test_value('Nearctic beak gape', nearctic_cities$gape_width_fdiv_normalised)

Warning: cannot compute exact p-value with ties

[1] "Nearctic beak gape mean 0.59 "

test_value('Nearctic Kipps distance', nearctic_cities$kipps_distance_fdiv_normalised)

Warning: cannot compute exact p-value with ties

[1] "Nearctic Kipps distance mean 0.58 "

Neartic Cluster 1`

nearactic_cluster1 = nearctic_cities %>% filter(cluster == 1)
plot_city_cluster(nearactic_cluster1, 'Neartic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster1.jpg'), width = get_image_width(nearactic_cluster1), height = get_image_height(nearactic_cluster1), units = "mm")

Neartic Cluster 2

nearactic_cluster2 = nearctic_cities %>% filter(cluster == 2)
plot_city_cluster(nearactic_cluster2, 'Neartic cluster 2')

Warning: cannot compute exact p-value with tiesWarning: cannot compute exact p-value with tiesWarning: cannot compute exact p-value with tiesWarning: cannot compute exact p-value with ties

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster2.jpg'), width = get_image_width(nearactic_cluster2), height = get_image_height(nearactic_cluster2), units = "mm")

Neartic Cluster 3

nearactic_cluster3 = nearctic_cities %>% filter(cluster == 3)
plot_city_cluster(nearactic_cluster3, 'Neartic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster3.jpg'), width = get_image_width(nearactic_cluster3), height = get_image_height(nearactic_cluster3), units = "mm")

Neartic Cluster 4

nearactic_cluster4 = nearctic_cities %>% filter(cluster == 4)
plot_city_cluster(nearactic_cluster4, 'Neartic cluster 4')

Warning: cannot compute exact p-value with ties

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster4.jpg'), width = get_image_width(nearactic_cluster4), height = get_image_height(nearactic_cluster4), units = "mm")

Neotropic

neotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Neotropic')
neotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
 [1] "Culiacán"                  "Guadalajara"               "Morelia"                   "Acapulco"                  "Querétaro"                
 [6] "Cuernavaca"                "Puebla"                    "Oaxaca"                    "Xalapa"                    "Veracruz"                 
[11] "Tuxtla Gutiérrez"          "Villahermosa"              "Guatemala City"            "San Salvador"              "San Pedro Sula"           
[16] "Mérida"                    "Tegucigalpa"               "Managua"                   "San José"                  "Cancún"                   
[21] "Guayaquil"                 "Chiclayo"                  "Panama City"               "Trujillo"                  "Quito"                    
[26] "Havana"                    "Cali"                      "Lima"                      "Pereira"                   "Miami"                    
[31] "Medellín"                  "Ibagué"                    "Cartagena"                 "Kingston"                  "Bogota"                   
[36] "Barranquilla"              "Bucaramanga"               "Cúcuta"                    "Maracaibo"                 "Arequipa"                 
[41] "Barquisimeto"              "Santo Domingo"             "Maracay"                   "El Alto [La Paz]"          "Caracas"                  
[46] "Cochabamba"                "Viña del Mar [Valparaíso]" "Río Piedras [San Juan]"    "Barcelona"                 "Concepción"               
[51] "Santiago"                  "Mendoza"                   "Salta"                     "Cordoba"                   "Asuncion"                 
[56] "Buenos Aires"              "La Plata"                  "Ciudad del Este"           "Montevideo"                "Mar del Plata"            
[61] "Porto Alegre"              "São Paulo"                 "Santos"                    "Sao Jose dos Campos"      

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

neotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% neotropic_cities_community_data$city_id)) 

Run 0 stress 0.134619 
Run 1 stress 0.134433 
... New best solution
... Procrustes: rmse 0.006824014  max resid 0.04572879 
Run 2 stress 0.1407144 
Run 3 stress 0.1348237 
... Procrustes: rmse 0.006646789  max resid 0.04604869 
Run 4 stress 0.1348046 
... Procrustes: rmse 0.006551306  max resid 0.04608394 
Run 5 stress 0.134433 
... Procrustes: rmse 0.00003293268  max resid 0.0001049247 
... Similar to previous best
Run 6 stress 0.134433 
... Procrustes: rmse 0.00005189183  max resid 0.0001643412 
... Similar to previous best
Run 7 stress 0.1366488 
Run 8 stress 0.1407144 
Run 9 stress 0.1348046 
... Procrustes: rmse 0.006551801  max resid 0.04608656 
Run 10 stress 0.134619 
... Procrustes: rmse 0.006841563  max resid 0.04571578 
Run 11 stress 0.1348237 
... Procrustes: rmse 0.006648345  max resid 0.04605669 
Run 12 stress 0.1344331 
... Procrustes: rmse 0.00006317981  max resid 0.0001865383 
... Similar to previous best
Run 13 stress 0.1405637 
Run 14 stress 0.1346406 
... Procrustes: rmse 0.007270798  max resid 0.04575267 
Run 15 stress 0.1405831 
Run 16 stress 0.1406945 
Run 17 stress 0.1346406 
... Procrustes: rmse 0.007285202  max resid 0.04573278 
Run 18 stress 0.1344509 
... Procrustes: rmse 0.001186371  max resid 0.006935107 
... Similar to previous best
Run 19 stress 0.1405634 
Run 20 stress 0.1344331 
... Procrustes: rmse 0.00007290489  max resid 0.0002305292 
... Similar to previous best
*** Best solution repeated 5 times

neotropic_cities_nmds

scree_plot(neotropic_cities_nmds)

neotropic_cities = cluster_cities(city_nmds = neotropic_cities_nmds, cities_community_data = neotropic_cities_community_data, centers = 5)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(neotropic_cities)

neotropic_biomes = resolve[resolve$REALM == 'Neotropic',c('REALM')]
 
ggplot() + 
  geom_sf(data = neotropic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = neotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_clusters.jpg'))

Saving 7.29 x 4.51 in image

test_value('Neotropic MNTD', neotropic_cities$mntd_normalised)

[1] "Neotropic MNTD mean 0.53 "

test_value('Neotropic beak gape', neotropic_cities$gape_width_fdiv_normalised)

[1] "Neotropic beak gape mean 0.44 "

test_value('Neotropic kipps distance', neotropic_cities$kipps_distance_fdiv_normalised)

[1] "Neotropic kipps distance mean 0.58 *"

Neotropic Cluster 1

neotropic_cluster1 = neotropic_cities %>% filter(cluster == 1)
plot_city_cluster(neotropic_cluster1, 'Neotropic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster1.jpg'), width = get_image_width(neotropic_cluster1), height = get_image_height(neotropic_cluster1), units = "mm")

Neotropic Cluster 2

neotropic_cluster2 = neotropic_cities %>% filter(cluster == 2)
plot_city_cluster(neotropic_cluster2, 'Neotropic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster2.jpg'), width = get_image_width(neotropic_cluster2), height = get_image_height(neotropic_cluster2), units = "mm")

Neotropic Cluster 3

neotropic_cluster3 = neotropic_cities %>% filter(cluster == 3)
plot_city_cluster(neotropic_cluster3, 'Neotropic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster3.jpg'), width = get_image_width(neotropic_cluster3), height = get_image_height(neotropic_cluster3), units = "mm")

Neotropic Cluster 4

neotropic_cluster4 = neotropic_cities %>% filter(cluster == 4)
plot_city_cluster(neotropic_cluster4, 'Neotropic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster4.jpg'), width = get_image_width(neotropic_cluster4), height = get_image_height(neotropic_cluster4), units = "mm")

Neotropic Cluster 5

neotropic_cluster5 = neotropic_cities %>% filter(cluster == 5)
plot_city_cluster(neotropic_cluster5, 'Neotropic cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster5.jpg'), width = get_image_width(neotropic_cluster5), height = get_image_height(neotropic_cluster5), units = "mm")

Palearctic

palearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Palearctic')
palearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
 [1] "Lisbon"                "Porto"                 "Marrakesh"             "Seville"               "Dublin"                "Málaga"               
 [7] "Madrid"                "Glasgow"               "Bilbao"                "Liverpool"             "Bristol"               "Manchester"           
[13] "Birmingham"            "Leeds"                 "Newcastle upon Tyne"   "Sheffield"             "Nottingham"            "Valencia"             
[19] "London"                "Toulouse"              "Paris"                 "Barcelona"             "Rotterdam [The Hague]" "Brussels"             
[25] "Amsterdam"             "Lyon"                  "Marseille"             "Dusseldorf"            "Nice"                  "Frankfurt am Main"    
[31] "Zurich"                "Oslo"                  "Stuttgart"             "Hamburg"               "Genoa"                 "Nuremberg"            
[37] "Copenhagen"            "Munich"                "Berlin"                "Dresden"               "Rome"                  "Prague"               
[43] "Stockholm"             "Poznan"                "Vienna"                "Wroclaw"               "Zagreb"                "Gdansk"               
[49] "Budapest"              "Krakow"                "Warsaw"                "Helsinki"              "Riga"                  "Belgrade"             
[55] "Lviv"                  "Sofia"                 "Thessaloniki"          "Saint Petersburg"      "Minsk"                 "Athens"               
[61] "Kyiv"                  "Istanbul"              "Odesa"                 "Samsun"                "Luxor"                 "Tel Aviv"             
[67] "Jerusalem"             "Tbilisi"               "Yerevan"               "Kuwait City"           "Doha"                  "Abu Dhabi"            
[73] "Dubai"                 "Bishkek"              

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

palearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% palearctic_cities_community_data$city_id)) 

Run 0 stress 0.04961857 
Run 1 stress 0.06522066 
Run 2 stress 0.06770805 
Run 3 stress 0.06505779 
Run 4 stress 0.08915955 
Run 5 stress 0.05395855 
Run 6 stress 0.04961847 
... New best solution
... Procrustes: rmse 0.0001353415  max resid 0.0006424586 
... Similar to previous best
Run 7 stress 0.07723898 
Run 8 stress 0.06982598 
Run 9 stress 0.05002297 
... Procrustes: rmse 0.008523462  max resid 0.02080155 
Run 10 stress 0.05001003 
... Procrustes: rmse 0.06742951  max resid 0.2228023 
Run 11 stress 0.05354024 
Run 12 stress 0.07364969 
Run 13 stress 0.04967779 
... Procrustes: rmse 0.02454439  max resid 0.1155996 
Run 14 stress 0.06518621 
Run 15 stress 0.0547908 
Run 16 stress 0.0675422 
Run 17 stress 0.06304729 
Run 18 stress 0.04968731 
... Procrustes: rmse 0.02483805  max resid 0.1165918 
Run 19 stress 0.06829409 
Run 20 stress 0.08070857 
*** Best solution repeated 1 times

palearctic_cities_nmds

scree_plot(palearctic_cities_nmds)

palearctic_cities = cluster_cities(city_nmds = palearctic_cities_nmds, cities_community_data = palearctic_cities_community_data, centers = 7)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(palearctic_cities)

palearctic_biomes = st_crop(resolve[resolve$REALM == 'Palearctic',c('REALM')], xmin = -30, ymin = 20, xmax = 80, ymax = 65)

although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning: attribute variables are assumed to be spatially constant throughout all geometries

 
ggplot() + 
  geom_sf(data = palearctic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = palearctic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_clusters.jpg'))

Saving 7.29 x 4.51 in image

test_value('Palearctic MNTD', palearctic_cities$mntd_normalised)

[1] "Palearctic MNTD mean 0.58 *"

test_value('Palearctic beak gape', palearctic_cities$gape_width_fdiv_normalised)

[1] "Palearctic beak gape mean 0.86 ***"

test_value('Palearctic kipps distance', palearctic_cities$kipps_distance_fdiv_normalised)

[1] "Palearctic kipps distance mean 0.91 ***"

Palearctic Cluster 1

palearctic_cluster1 = palearctic_cities %>% filter(cluster == 1)
plot_city_cluster(palearctic_cluster1, 'Palearctic cluster 1')

Warning: cannot compute exact p-value with ties

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster1.jpg'), width = get_image_width(palearctic_cluster1), height = get_image_height(palearctic_cluster1), units = "mm")

Palearctic Cluster 2

palearctic_cluster2 = palearctic_cities %>% filter(cluster == 2)
plot_city_cluster(palearctic_cluster2, 'Palearctic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster2.jpg'), width = get_image_width(palearctic_cluster2), height = get_image_height(palearctic_cluster2), units = "mm")

Palearctic Cluster 3

palearctic_cluster3 = palearctic_cities %>% filter(cluster == 3)
plot_city_cluster(palearctic_cluster3, 'Palearctic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster3.jpg'), width = get_image_width(palearctic_cluster3), height = get_image_height(palearctic_cluster3), units = "mm")

Palearctic Cluster 4

palearctic_cluster4 = palearctic_cities %>% filter(cluster == 4)
plot_city_cluster(palearctic_cluster4, 'Palearctic cluster 4')

Warning: cannot compute exact p-value with tiesWarning: cannot compute exact p-value with ties

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster4.jpg'), width = get_image_width(palearctic_cluster4), height = get_image_height(palearctic_cluster4), units = "mm")

Palearctic Cluster 5

palearctic_cluster5 = palearctic_cities %>% filter(cluster == 5)
plot_city_cluster(palearctic_cluster5, 'Palearctic cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster5.jpg'), width = get_image_width(palearctic_cluster5), height = get_image_height(palearctic_cluster5), units = "mm")

Palearctic Cluster 6

palearctic_cluster6 = palearctic_cities %>% filter(cluster == 6)
plot_city_cluster(palearctic_cluster6, 'Palearctic cluster 6')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster6.jpg'), width = get_image_width(palearctic_cluster6), height = get_image_height(palearctic_cluster6), units = "mm")

Palearctic Cluster 7

palearctic_cluster7 = palearctic_cities %>% filter(cluster == 7)
plot_city_cluster(palearctic_cluster7, 'Palearctic cluster 7')

Warning: cannot compute exact p-value with ties

ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster7.jpg'), width = get_image_width(palearctic_cluster7), height = get_image_height(palearctic_cluster7), units = "mm")

Afrotropic

afrotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Afrotropic')
afrotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
[1] "Cape Town"    "Johannesburg" "Pretoria"     "Kigali"       "Kampala"      "Arusha"       "Nairobi"      "Addis Ababa"  "Antananarivo"

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

afrotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% afrotropic_cities_community_data$city_id)) 

Run 0 stress 0.00009014786 
Run 1 stress 0.0000894423 
... New best solution
... Procrustes: rmse 0.0002302813  max resid 0.0003204519 
... Similar to previous best
Run 2 stress 0.0002552434 
... Procrustes: rmse 0.001471005  max resid 0.001947622 
... Similar to previous best
Run 3 stress 0.00009602873 
... Procrustes: rmse 0.0003134658  max resid 0.0004532654 
... Similar to previous best
Run 4 stress 0.0004206102 
... Procrustes: rmse 0.002382671  max resid 0.003350133 
... Similar to previous best
Run 5 stress 0.001267824 
Run 6 stress 0.00009984846 
... Procrustes: rmse 0.000573798  max resid 0.0007724523 
... Similar to previous best
Run 7 stress 0.00009404138 
... Procrustes: rmse 0.0002198999  max resid 0.0002986639 
... Similar to previous best
Run 8 stress 0.00006937968 
... New best solution
... Procrustes: rmse 0.0001548829  max resid 0.0003491895 
... Similar to previous best
Run 9 stress 0.00009515083 
... Procrustes: rmse 0.0002518215  max resid 0.0003915592 
... Similar to previous best
Run 10 stress 0.0002810887 
... Procrustes: rmse 0.001624066  max resid 0.002239305 
... Similar to previous best
Run 11 stress 0.3209238 
Run 12 stress 0.0001607375 
... Procrustes: rmse 0.0009389978  max resid 0.001284128 
... Similar to previous best
Run 13 stress 0.00009851387 
... Procrustes: rmse 0.0002528255  max resid 0.0003960776 
... Similar to previous best
Run 14 stress 0.00009299962 
... Procrustes: rmse 0.0002384291  max resid 0.0003789721 
... Similar to previous best
Run 15 stress 0.0001360544 
... Procrustes: rmse 0.0007908028  max resid 0.001108713 
... Similar to previous best
Run 16 stress 0.00009817612 
... Procrustes: rmse 0.0002561838  max resid 0.0003936572 
... Similar to previous best
Run 17 stress 0.002120919 
Run 18 stress 0.00009438058 
... Procrustes: rmse 0.000245643  max resid 0.0003771271 
... Similar to previous best
Run 19 stress 0.00009470345 
... Procrustes: rmse 0.0002746505  max resid 0.0004407572 
... Similar to previous best
Run 20 stress 0.00008594013 
... Procrustes: rmse 0.0002324527  max resid 0.0003605282 
... Similar to previous best
*** Best solution repeated 11 times

Warning: stress is (nearly) zero: you may have insufficient data

afrotropic_cities_nmds

scree_plot(afrotropic_cities_nmds)

afrotropic_cities = cluster_cities(city_nmds = afrotropic_cities_nmds, cities_community_data = afrotropic_cities_community_data, centers = 2)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(afrotropic_cities)

afrotropic_biomes = resolve[resolve$REALM == 'Afrotropic',c('REALM')]
 
ggplot() + 
  geom_sf(data = afrotropic_biomes, aes(geometry = geometry)) + 
  geom_sf(data = afrotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_clusters.jpg'))

Saving 7.29 x 4.51 in image

test_value('Afrotropic MNTD', afrotropic_cities$mntd_normalised)

[1] "Afrotropic MNTD mean 0.14 *"

test_value('Afrotropic beak gape', afrotropic_cities$gape_width_fdiv_normalised)

[1] "Afrotropic beak gape mean 0.41 "

test_value('Afrotropic kipps distance', afrotropic_cities$kipps_distance_fdiv_normalised)

[1] "Afrotropic kipps distance mean 0.42 "

Afrotropic Cluster 1

afrotropic_cluster1 = afrotropic_cities %>% filter(cluster == 1)
plot_city_cluster(afrotropic_cluster1, 'Afrotropic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster1.jpg'), width = get_image_width(afrotropic_cluster1), height = get_image_height(afrotropic_cluster1), units = "mm")

Afrotropic Cluster 2

afrotropic_cluster2 = afrotropic_cities %>% filter(cluster == 2)
plot_city_cluster(afrotropic_cluster2, 'Afrotropic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster2.jpg'), width = get_image_width(afrotropic_cluster2), height = get_image_height(afrotropic_cluster2), units = "mm")

Indomalayan

indomalayan_cities_community_data = community_data_metrics %>% filter(core_realm == 'Indomalayan')
indomalayan_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()

$city_name
  [1] "Srinagar"            "Jamnagar"            "Jammu"               "Rajkot"              "Bikaner"             "Jodhpur"             "Jalandhar"          
  [8] "Ahmedabad"           "Bhavnagar"           "Ludhiana"            "Anand"               "Udaipur"             "Surat"               "Vadodara"           
 [15] "Ajmer"               "Chandigarh"          "Vasai-Virar"         "Mumbai"              "Jaipur"              "Delhi [New Delhi]"   "Nashik"             
 [22] "Dehradun"            "Kota"                "Pune"                "Haridwar"            "Dhule"               "Ujjain"              "Indore"             
 [29] "Ahmadnagar"          "Kolhapur"            "Jalgaon"             "Agra"                "Aurangabad"          "Sangli"              "Belagavi"           
 [36] "Gwalior"             "Budaun"              "Bareilly"            "Dharwad"             "Bhopal"              "Bhind"               "Mangaluru"          
 [43] "Solapur"             "Vijayapura"          "Akola"               "Latur"               "Kannur"              "Davanagere"          "Thalassery"         
 [50] "Amravati"            "Kalaburagi"          "Kozhikode"           "Guruvayur"           "Malappuram"          "Lucknow"             "Thrissur"           
 [57] "Mysuru"              "Kochi"               "Alappuzha"           "Nagpur"              "Kollam"              "Jabalpur"            "Ettumanoor"         
 [64] "Hyderabad"           "Coimbatore"          "Bengaluru"           "Thiruvananthapuram"  "Tiruppur"            "Faizabad"            "Erode"              
 [71] "Prayagraj"           "Pratapgarh"          "Salem"               "Dindigul"            "Madurai"             "Tiruchirappalli"     "Durg"               
 [78] "Vellore"             "Tirupati"            "Raipur"              "Bilaspur"            "Vijayawada"          "Puducherry"          "Chennai"            
 [85] "Kathmandu"           "Colombo"             "Rajamahendravaram"   "Patna"               "Kandy"               "Bihar Sharif"        "Visakhapatnam"      
 [92] "Ranchi"              "Brahmapur"           "Jamshedpur"          "Darjeeling"          "Siliguri"            "Cuttack"             "Bhubaneshwar"       
 [99] "Jalpaiguri"          "Berhampore"          "Kolkata"             "Krishnanagar"        "Guwahati [Dispur]"   "Agartala"            "Silchar"            
[106] "Dimapur"             "Bangkok"             "George Town"         "Kuala Lumpur"        "Phnom Penh"          "Singapore"           "Hong Kong"          
[113] "Sha Tin"             "Hsinchu"             "Taichung"            "New Taipei [Taipei]" "Tainan"              "Denpasar"            "Kaohsiung"          
[120] "Kota Kinabalu"      

attr(,"na.action")
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
  1  56  81  87  90  92  94  96  98 100 102 103 104 106 112 115 116 117 118 119 121 215 
attr(,"class")
[1] "omit"

indomalayan_cities_nmds = community_nmds(communities %>% filter(city_id %in% indomalayan_cities_community_data$city_id)) 

Run 0 stress 0.1190668 
Run 1 stress 0.1224401 
Run 2 stress 0.1394961 
Run 3 stress 0.1161598 
... New best solution
... Procrustes: rmse 0.02798297  max resid 0.2346031 
Run 4 stress 0.1540252 
Run 5 stress 0.1175538 
Run 6 stress 0.1167657 
Run 7 stress 0.1384412 
Run 8 stress 0.1241275 
Run 9 stress 0.1163498 
... Procrustes: rmse 0.025336  max resid 0.2481842 
Run 10 stress 0.1153501 
... New best solution
... Procrustes: rmse 0.008287738  max resid 0.08217213 
Run 11 stress 0.1199235 
Run 12 stress 0.1172201 
Run 13 stress 0.151756 
Run 14 stress 0.1246307 
Run 15 stress 0.134023 
Run 16 stress 0.1580175 
Run 17 stress 0.1229198 
Run 18 stress 0.1502633 
Run 19 stress 0.1189366 
Run 20 stress 0.117072 
Run 21 stress 0.1458634 
Run 22 stress 0.137476 
Run 23 stress 0.1186188 
Run 24 stress 0.1637221 
Run 25 stress 0.1191928 
Run 26 stress 0.1164559 
Run 27 stress 0.142366 
Run 28 stress 0.1379337 
Run 29 stress 0.1191853 
Run 30 stress 0.1266232 
*** Best solution was not repeated -- monoMDS stopping criteria:
    29: stress ratio > sratmax
     1: scale factor of the gradient < sfgrmin

indomalayan_cities_nmds

scree_plot(indomalayan_cities_nmds)

indomalayan_cities = cluster_cities(city_nmds = indomalayan_cities_nmds, cities_community_data = indomalayan_cities_community_data, centers = 5)

Joining with `by = join_by(city_id)`

plot_nmds_clusters(indomalayan_cities)

indomalayan_biomes = resolve[resolve$REALM == 'Indomalayan',c('REALM')]
 
ggplot() + 
  geom_sf(data = indomalayan_biomes, aes(geometry = geometry)) + 
  geom_sf(data = indomalayan_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_clusters.jpg'))

Saving 7.29 x 4.51 in image

test_value('Indomalayan MNTD', indomalayan_cities$mntd_normalised)

[1] "Indomalayan MNTD mean 0.44 ***"

test_value('Indomalayan beak gape', indomalayan_cities$gape_width_fdiv_normalised)

[1] "Indomalayan beak gape mean 0.52 "

test_value('Indomalayan kipps distance', indomalayan_cities$kipps_distance_fdiv_normalised)

[1] "Indomalayan kipps distance mean 0.82 ***"

Indomalayan Cluster 1

indomalayan_cluster1 = indomalayan_cities %>% filter(cluster == 1)
plot_city_cluster(indomalayan_cluster1, 'Indomalayan cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster1.jpg'), width = get_image_width(indomalayan_cluster1), height = get_image_height(indomalayan_cluster1), units = "mm")

Indomalayan Cluster 2

indomalayan_cluster2 = indomalayan_cities %>% filter(cluster == 2)
plot_city_cluster(indomalayan_cluster2, 'Indomalayan cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster2.jpg'), width = get_image_width(indomalayan_cluster2), height = get_image_height(indomalayan_cluster2), units = "mm")

Indomalayan Cluster 3

indomalayan_cluster3 = indomalayan_cities %>% filter(cluster == 3)
plot_city_cluster(indomalayan_cluster3, 'Indomalayan cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster3.jpg'), width = get_image_width(indomalayan_cluster3), height = get_image_height(indomalayan_cluster3), units = "mm")

Indomalayan Cluster 4

indomalayan_cluster4 = indomalayan_cities %>% filter(cluster == 4)
plot_city_cluster(indomalayan_cluster4, 'Indomalayan cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster4.jpg'), width = get_image_width(indomalayan_cluster4), height = get_image_height(indomalayan_cluster4), units = "mm")

Indomalayan Cluster 5

indomalayan_cluster5 = indomalayan_cities %>% filter(cluster == 5)
plot_city_cluster(indomalayan_cluster5, 'Indomalayan cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster5.jpg'), width = get_image_width(indomalayan_cluster5), height = get_image_height(indomalayan_cluster5), units = "mm")